Journal  of  Power  Sources  207  (2012)  150-159 


ELSEVIER 


Contents  lists  available  at  SciVerse  ScienceDirect 

Journal  of  Power  Sources 

journal  homepage:  www.elsevier.com/locate/jpowsour 


in 

SllbudtS 


A  second  nearest-neighbor  embedded  atom  method  interatomic  potential  for 
Li-Si  alloys 

Zhiwei  Cuia,  Feng  Gaoa,  Zhihua  Cuib,  Jianmin  Qua* 

a  Department  of  Civil  and  Environmental  Engineering,  Northwestern  University,  Evanston,  ll  60208,  USA 

b  Complex  System  and  Computational  Intelligence  Laboratory,  Taiyuan  University  of  Science  and  Technology,  Taiyuan,  Shanxi  Province  030024,  PR  China 


ARTICLE  INFO 


ABSTRACT 


Article  history: 

Received  21  October  2011 

Received  in  revised  form  28  January  2012 

Accepted  29  January  2012 

Available  online  6  February  2012 


Keywords: 

Li-ion  battery 
Li-Si  alloy 

Modified  embedded  atom  method 
Particle  swarm  optimization 
Disordered-ordered  transition 


A  second  nearest-neighbor  modified  embedded  atom  method  (2NN  MEAM)  interatomic  potential  for 
lithium-silicon  (Li-Si)  alloys  is  developed  by  using  the  particle  swarm  optimization  (PSO)  method  in 
conjunction  with  ab  initio  calculations.  It  is  shown  that  the  new  interatomic  potential  is  capable  of  sim¬ 
ulating  the  transition  from  disordered  to  ordered  states  of  Li-Si  crystalline  structures,  an  indication  of 
the  stability  and  robustness  of  the  interatomic  potential  at  finite  temperature.  Examples  are  given  to 
demonstrate  that  the  new  interatomic  potential  is  also  capable  of  predicting  the  material  properties  of 
both  crystalline  and  amorphous  Li-Si  alloys,  including  the  elastic  modulus,  compositional  expansion, 
diffusivity  of  Li  in  Li-Si  alloys,  plastic  yield  strength,  etc. 
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1.  Introduction 

Lithium  (Li)  ion  batteries  have  gained  tremendous  attentions 
recently  in  high  energy  capacity,  high  operating  voltage  and  low 
self-discharge  energy  storage  devices  [1].  For  higher  capacity  Li- 
ion  batteries,  silicon  (Si)  is  being  considered  as  a  promising  anode 
material  owing  to  its  highest  known  theoretical  charge  capacity 
(4200  mAh  g-1)  [2,3].  However,  the  volume  expansion  of  the  Si 
anode  can  be  as  much  as  400%  after  fully  charged,  which  often 
leads  fracture  failure  resulting  in  reduced  battery  life  [4-6].  Such 
mechanical  failure  has  been  one  of  the  major  roadblocks  that 
prevent  the  technology  breakthroughs  for  large  scale  commercial¬ 
ization  of  Si-based  anode  Li-ion  batteries. 

Extensive  numerical  simulations  [7-9]  and  microscopic  ana¬ 
lyzes  [4,6,10-19]  have  confirmed  that  the  insertion  of  Li  ions  into 
Si  forms  Li-Si  alloys.  Depending  on  the  amount  of  Li,  Li*Si  alloys  of 
various  compositions  might  be  formed  [10,11].  Volumetric  expan¬ 
sion  due  to  Li  insertion  has  also  been  studied  extensively  both 
numerically  [7,20]  and  experimentally  [8,21-23].  Based  on  known 
compositions  and  their  volume  change  with  respect  to  Si  lattice, 
many  researchers  [24-32]  have  computed  the  stresses  generated 
during  Li  insertion. 

Although  many  attempts  have  been  made  and  several  models 
have  been  developed  for  understanding  the  precise  mechanisms  of 
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Li  insertion/extraction,  these  models  are  still  far  from  being  quanti¬ 
tative  and  are  not  ready  for  design  applications.  For  instance,  some 
existing  models  predict  the  stress  generated  by  Li  insertion,  but  do 
not  account  for  the  effect  of  stress  on  the  diffusion  of  Li.  Without 
such  two-way  coupling,  these  exiting  models  are  unable  to  simu¬ 
late  correctly  the  insertion  of  Li  in  the  Si,  thus  cannot  calculate  the 
Li  concentration  accurately.  Without  such  information,  one  cannot 
calculate  the  stresses  correctly,  nor  can  one  predict  the  fracture  of 
the  Si  anode.  This  is  because  Li  concentration  is  non-uniform,  and 
highly  dependent  on  the  diffusion  of  Li  within  the  Si  anode,  which  is 
strongly  influenced  by  stress  distribution  [21,33].  Although  some 
recent  studies  have  attempted  to  investigate  the  effect  of  stress 
on  Li  diffusion  [32-36],  models  used  in  these  studies  either  did 
not  account  for  the  electrochemical  driving  forces  for  Li  transport 
[37,38  ],  or  ignored  the  effect  of  shear  stresses,  or  was  based  on  small 
strain  deformation. 

More  importantly,  recent  TEM  in  situ  observations  [39-41  ]  show 
that  during  the  first  charge  of  a  Si  nanoparticle  or  nanowire,  a  sharp 
crystal/amorphous  interface  is  formed  that  moves  towards  the  cen¬ 
ter  as  the  charge  progresses.  At  some  point,  fracture  occurs  on  the 
particle  surface  due  to  tensile  hoop  stress.  However,  none  of  the 
existing  theories  seems  to  be  able  to  predict  such  failure  process. 

A  key  drawback  of  all  the  existing  models  is  the  lack  of  input 
from  material  composition  and  microstructure.  In  fact,  all  the  avail¬ 
able  models  are  based  on  continuum  mechanics.  Yet,  Li  insertion 
is  inherently  atomistic,  because  the  insertion  of  Li  atoms  into  Si 
causes  significant  changes  in  the  Si  bond  structure  [42].  What 
this  means  is  that  care  must  be  taken  in  developing  continuum 
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models  in  order  to  account  for  such  microstructural  changes  during 
insertion.  Otherwise,  continuum  models  would  be  unable  to  simu¬ 
late  certain  physical  phenomena  such  as  phase  transformation  and 
fracture  as  observed  in  Li  insertion. 

Another  major  shortcoming  of  the  existing  continuum  models 
is  the  lack  of  appropriate  values  for  many  of  the  physical  proper¬ 
ties  in  the  model,  such  as  the  plastic  yield  strength  of  LixSi  as  a 
function  of  Li  concentration,  and  the  diffusivity  of  Li  in  LixSi  and  its 
dependence  on  the  stress.  In  theory,  these  properties  can  be  mea¬ 
sured  experimentally.  In  reality,  such  measurements  are  extremely 
cumbersome,  and  difficult  to  perform. 

With  recent  advances  in  computational  algorithms  and  com¬ 
puter  hardware,  atomistic  level  computations,  such  as  molecular 
dynamic  (MD)  simulations  have  become  an  effective  and  efficient 
tool  for  simulating  material  behavior  and  for  computing  material 
properties.  If  done  correctly,  MD  simulations  can  provide  insightful 
information  on  the  microstructural  changes  during  deformation, 
and  enhance  our  understanding  of  the  failure  mechanisms.  We 
believe  that  many  of  the  existing  continuum  models  on  Li  inser¬ 
tion  will  benefit  greatly  from  comparing  with  atomistic  level 
simulations. 

However,  to  conduct  MD  simulations,  the  first  thing  needed  is 
a  proper  interatomic  potential.  At  the  present  time,  interatomic 
potentials  for  Li-Si  alloys  are  not  available.  This  is  the  major  reason 
that  no  MD  simulation  results  have  been  published  in  the  open 
literature  so  far.  In  fact,  Haftbaradaran  et  al.  [34]  had  developed 
a  continuum  model  for  Li  diffusion  in  LixSi  alloys.  To  verify  their 
continuum  model,  they  have  conducted  MD  simulations.  However, 
because  of  the  lack  of  interatomic  potential  for  Li-Si  alloys,  their 
MD  simulations  were  performed  on  the  diffusion  of  hydrogen  in 
single  crystal  fee  nickel. 

Clearly,  there  is  a  critical  need  for  developing  an  accurate  and 
robust  interatomic  potential  for  conducting  MD  simulations  of  Li 
insertion  into  Si  anode.  The  potential  should  be  able  to  describe  all 
the  possible  atomic  interactions  in  all  the  possible  crystal  struc¬ 
tures  of  LixSi,  as  well  as  the  amorphous  LixSi  of  any  0<x<4.4. 
It  should  also  provide  stable  solution  at  both  room  and  elevated 
temperatures.  Furthermore,  the  potential  must  be  able  to  calcu¬ 
late  accurately  the  basic  properties  of  LixSi  including  the  elastic 
constants,  the  lattice  constants,  the  cohesive  energies,  etc. 

Such  an  interatomic  potential  for  LixSi  will  be  a  powerful  tool.  It 
enables  us  to  validate  the  various  continuum  models  on  Li  insertion. 
It  can  be  used  to  estimate  some  of  the  material  properties  that  are 
difficult  to  obtain  experimentally.  It  gives  us  a  mean  to  study  the 
microstructural  mechanisms  of  fracture  and  failure. 

One  of  the  most  widely  used  interatomic  potentials  for  metals  is 
the  modified  embedded  atom  method  (MEAM)  potential  first  pro¬ 
posed  by  Baskes  [43  ].  It  has  been  used  for  a  variety  of  fee  metals  with 
great  success.  However,  as  pointed  out  recently  by  Lee  et  al.  [44], 
the  MEAM  potential  has  some  drawbacks  when  used  for  bcc  crystal 
structures.  For  instance,  the  (111)  surface  energy  was  found  to  be 
smaller  than  that  of  (1  0  0)  surface  by  the  MEAM  potential,  which  is 
contrary  to  experimental  results  or  ab  initio  calculations.  Also,  the 
MEAM  potential  predicts  that  the  most  stable  structure  for  Li  is  not 
bcc,  which  is  obviously  incorrect. 

To  better  describe  the  bcc  structure,  Lee  et  al.  [44]  proposed  the 
second  nearest  neighbor  (2NN)  MEAM  interatomic  potential,  and 
demonstrated  that  it  is  better  suited  for  the  bcc  structure.  Since 
then,  2NN  MEAM  potentials  have  been  developed  for  several  ele¬ 
ments  and  their  alloys  including  silicon  [45]  and  magnesium  [46]. 

This  paper  describes  the  development  of  a  2NN  MEAM  poten¬ 
tial  for  the  LixSi  alloys.  Our  general  approach  is  to  use  data  from 
existing  experimental  measurements  and  ab  initio  computations. 
Based  on  these  data,  the  parameters  in  the  2NN  MEAM  potential 
will  be  optimized  by  the  particle  swarm  optimization  (PSO)  [47-49] 
technique. 


Before  proceeding,  however,  we  note  that  LixSi  alloys  are  some¬ 
what  different  from  typical  metal  alloys  in  that  they  exhibit  some 
behavior  of  ionic  solids.  For  example,  a  number  of  investigators 
[50,51]  have  used  first  principles  calculations  to  study  the  charge 
transfer  in  LixSi  crystal  structures.  Their  results  reveal  that  Li  atoms 
donate  a  certain  amount  of  electrons  to  the  Si  atoms  in  all  crystalline 
Li-Si  phases.  However,  the  charge  carried  by  the  Si  atoms  seems  to 
differ  depending  on  the  Li  concentration,  and  the  crystalline  struc¬ 
ture.  Furthermore,  in  a  given  Li-Si  crystalline  structure,  the  charge 
distribution  among  Si  atoms  may  not  be  the  same.  This  poses  a  fun¬ 
damental  question  as  to  whether  an  interatomic  potential  exists 
that  is  universal  to  all  the  different  compositions  and  structures  of 
LixSi.  Although  the  answer  to  this  question  is  not  clear  at  this  point, 
our  preliminary  MD  results,  as  discussed  in  this  paper,  do  show 
good  agreement  with  the  experimental  data  or  with  the  ab  initio 
computations,  indicating  that  charge  transfer  effect  in  LixSi  may  be 
neglected,  and  an  accurate  and  robust  interatomic  potential  might 
exist  that  can  accurately  predict  the  behavior  of  LixSi  alloys  under 
a  wide  range  of  conditions. 


2.  Methodology 

2.1.  2NN  MEAM  potential 


In  a  2NN  MEAM  formulation,  the  total  energy  of  a  system  is 
written  as  [45,52] 


rcpi) + \y^<Pij(Rn) 


K  *  0 


(1) 


where  F2(pj)  is  the  embedding  function,  p2  is  the  background  elec¬ 
tron  density  at  the  site  where  atom  i  occupies,  and  (pij{Ry)  is  the  pair 
interaction  between  atoms  i  and  j  at  a  distance  Ry.  The  background 
electron  density  p  is  composed  of  several  partial  electron  density 
terms.  Each  partial  electron  density  is  a  function  of  atomic  config¬ 
uration  and  atomic  electron  density.  The  atomic  electron  density  is 
given  by 


pfh\Rij)  =  Poexp 


(2) 


where  p0  is  a  scaling  factor,  are  adjustable  parameters  and  re  is 
the  nearest  neighbor  distance  in  the  reference  structure  after  equi¬ 
librium.  In  the  2NN  MEAM,  the  energy  per  atom  in  the  reference 
structure  is  estimated  from  a  universal  equation  of  state  [53].  Once 
Fj(pj)  is  known,  the  pair  potential  (Pij(Ry)  can  be  evaluated.  More 
details  of  the  MEAM  formulation  can  be  found  in  [45,52]. 

For  a  binary  intermetallic  compound,  the  only  unknown  func¬ 
tion  in  the  potential  is  the  pair  potential  function  between  the 
different  types  of  atoms.  To  obtain  the  pair  potential  function, 
the  atoms  that  have  the  same  type  of  atoms  as  second  nearest- 
neighbors  are  considered  as  a  reference  structure.  For  the  LixSi 
alloy  considered  here,  we  used  a  Li3 Si-type  Ll2  ordered  structure 
as  the  reference  structure.  Thus,  the  total  energy  per  atom  (for  3/4 
Li  atom  + 1  /4  Si  atom)  is  given  by, 


£Li3Si(R)=  -Ahi{Psi)+^u{Pu)+f 


2^SiLi(^)  +  2^LiLiW 


2 


1  3 

^5Si^SiSi(a^)  +  ^U<Puu{a^) 


(3) 


where  Z\  and  Z2  are  the  numbers  of  first  and  second  NN,  respec¬ 
tively,  in  the  Ll2  Li3Si  structure,  SLi  and  Ssi  are  the  screening 
function  for  the  2NN  interactions  between  Li,  and  Si  atoms,  respec¬ 
tively,  and  a  is  the  ratio  between  the  second  and  first  NN  distances 
in  the  reference  structure.  The  pair  potential  function  between  Li 
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and  Si  can  then  be  evaluated  by 

(PUSiW  =  ^£Li3Si(R)  -  ^Fsi(Psi)  -  \MPli)  -  ^LiLi(R) 

-  ^SSi^siSi(Q^)  -  ^SLi^LiLj(aR).  (4) 

The  only  unknown  term  on  the  right  hand  side  of  (4)  is  the  total 
energy  per  atom,  E^si.  To  obtain  E“^s[,  the  universal  equation  of 
state  can  be  used.  For  the  Ll2  structure  of  Li3Si, 

E“i3Si(R)  =  -Ec[l  +a*  +d(a*)3]e-a*,  (5) 

where 


In  the  above,  re  is  the  equilibrium  NN  distance,  d  is  an  adjustable 
parameter,  B  is  the  bulk  modulus,  and  Q  is  the  equilibrium 
atomic  volume.  Once  these  parameters  are  determined  for  the  Ll2 
structure,  the  corresponding  pair  potential  function  (pLis[  can  be 
evaluated  from  (4). 

In  addition  to  these  4  parameters  in  the  universal  equation  of 
state,  there  are  9  other  parameters  in  the  potential  that  need  to  be 
determined.  One  is  the  atomic  electron  density  scaling  factor  p0.  It 
has  been  found  that  p0  has  little  impact  on  pure  element,  but  plays 
a  significant  role  in  alloys  [44,45].  Furthermore,  Cmax  and  Cmin  are 
also  crucial  in  determining  the  interaction  range  of  the  alloys.  For 
Li-Si  binary  alloy  system,  there  are  4  combinations  of  both  Cmax  and 
Cmin.  i-e-  (Li— Si— Li,  Li— Li— Si,  Si— Si— Li,  Si— Li— Si)  for  a  total  of  8  param¬ 
eters.  Altogether,  there  are  13  parameters  that  can  be  adjusted  in 
order  to  fit  the  2NN  MEAM  potential  for  a  binary  alloy.  They  are  ECt 
re,  a,  d,  po,  and  4  combinations  of  both  Cmax  and  Cmin.  For  a  given 
set  of  these  13  parameters,  a  number  of  physical  properties  of  the 
five  known  crystalline  structures  of  Li*Si  alloys  (x  =  l,  12/7,  13/4, 
15/4, 22/5),  such  as  the  lattice  constants,  cohesive  energies,  elastic 
constants,  etc.,  can  be  predicted  by  MD  simulations.  By  fitting  these 
MD  predictions  to  their  corresponding  values  from  either  exper¬ 
imental  measurements  or  ab  initio  simulations,  the  13  adjustable 
parameters  in  the  2NN  MEAM  potential  can  be  optimized  to  yield 
an  accurate  and  hopefully  robust  interatomic  potential  for  the  Li-Si 
alloys.  In  the  next  section,  we  will  describe  how  such  optimization 
can  be  accomplished  by  the  PSO  method. 

To  obtain  the  complete  2NN  MEAM  potential  for  the  Li*Si  alloys, 
the  corresponding  potentials  for  Li  and  Si  are  also  needed.  The  2NN 
MEAM  potential  for  pure  Si  has  been  obtained  by  Lee  [45].  We  have 
modified  Lee’s  potential  slightly  so  it  is  better  suited  for  develop¬ 
ing  the  2NN  MEAM  potential  for  the  Li-Si  system,  see  Appendix  A. 
The  2NN  MEAM  potential  for  Li  was  also  developed  earlier  by  the 
present  authors  [52]. 

2.2.  Optimization  procedures 

In  this  section,  we  describe  how  the  particle  swarm  optimization 
(PSO)  method  [47-49]  can  be  used  to  determine  simultaneously 
the  13  adjustable  parameters  in  the  2NN  MEAM  potential.  The  PSO 
method  is  a  novel  population-based  optimization  technique.  It  was 
developed  initially  to  simulate  animal  social  behaviors,  e.g.,  birds 
flocking,  fish  schooling  and  insects  herding.  Due  to  its  simple  algo¬ 
rithm  and  fast  convergence,  the  PSO  has  been  widely  adopted  in 
many  areas  such  as  power  system  [54],  structural  damage  identi¬ 
fication  [55],  nonlinear  system  identification  [56],  ice-storage  air 
conditioning  system  [57],  etc. 

The  PSO  method  starts  with  a  group  (called  swarm)  of  candidate 
solutions  (called  particles).  These  particles  move  around  within  the 
search  space  to  seek  food  (called  optimum).  Let  77  (j)  be  an  objective 
function  that  transforms  a  particle  to  a  unique  real  number.  The  goal 


Table  1 

Optimized  parameters  in  the  2NN  MEAM  potential  for  Li-Si  alloys. 


Parameters 

Ec  (eVatorrr1) 

2.45 

re  (A) 

2.75 

a 

4.10 

d 

0.10 

Po  (Po/Po  ) 

3.0 

Cmax  (Li— Li— Si) 

2.81 

Cmax  (Si— Si— Li) 

2.20 

Cmax  (Li-Si— Li) 

2.40 

Cmax  (Li-Si-Si) 

2.40 

Cmin  (Li— Li— Si) 

0.55 

Cmin  (Si  Si— Li) 

0.35 

Cmin  (Li  Si— Li) 

0.45 

Cmin  (Li-Si-Si) 

0.45 

is  to  find  a  particular  particle  j  such  that  77(j)  <  /T(/<)  for  all  particle 
k  within  the  search  space. 

Let  x  =  (xi ,  X\ , . . .,  xm)  represent  a  set  of  candidate  parameters  in 
the  2NN  MEAM  for  a  specific  material.  Thus  the  objective  function 
for  the  optimization  can  be  defined  as, 


N 

77(x)  =  jy, 
1=1 


-fiW 
.  y\ 


(7) 


where  /j(x)  is  a  physical  property  calculated  by  MD  simulations 
using  the  2NN  MEAM  potential  with  the  parameter  set  x.  Such 
physical  properties  can  be,  for  example,  lattice  constants  or  elastic 
constants.  The  y2-  is  the  corresponding  value  of  the  same  physical 
property  from  ab  initio  calculations,  or  experimental  measure¬ 
ments.  The  weighting  factor  <r2  is  a  positive  number  selected  based 
on  the  importance  of  this  particular  physical  property.  In  this  study, 
structural  properties  are  considered  much  more  important  than 
auxiliary  structures,  so  they  are  given  much  larger  weights.  The 
objective  here  is  to  find  a  particular  x0  so  that  77(x0)  <  77(x)  for  all 
x  in  the  search  space. 

To  start  the  optimization,  each  parameter  is  assigned  a  random 
value  as  the  initial  input.  Physical  properties  of  the  systems  are 
then  computed  based  on  the  2NN  MEAM  potential  with  this  initial 
set  of  parameters.  This  gives  the  first  iteration  of  //(x)  in  Eq.  (7). 
The  next  set  of  candidate  parameters  is  then  selected  using  the 
PSO  algorithm.  A  small  enough  “time  step”  is  used  in  the  PSO  to 
ensure  that  each  partial  electron  density  term  has  the  same  order  of 
magnitude.  The  above  procedures  are  repeated  until  the  objective 
function  77(x)  has  been  minimized  to  a  satisfactory  level.  All  the 
MD  simulations  in  this  study  are  conducted  by  using  the  LAMMPS 
software  code  [58,59]. 

As  seen  from  Eq.  (7),  a  database  of  physical  properties  y2  is 
needed  in  order  to  optimize  the  parameters  in  the  2NN  MEAM 
potential.  In  this  work,  the  physical  properties  used  include  the 
lattice  constants,  cohesive  energies  and  elastic  constants  for  5 
LixSi  alloy  structures  (x  =  l,  12/7,  13/4,  15/4,  22/5,  corresponding 
to  Lii Sii ,  Lii2Si7,  Li^SU,  LiisSU,  and  L22Sis,  respectively),  as  well 
as  the  binding  energies  of  five  auxiliary  structures,  namely,  four 
crystals  (Bl,  B2,  B3,  LI 2)  and  one  molecular  structures  (LiSi).  These 
properties  are  either  calculated  by  ab  initio  simulations  carried  in 
this  work  or  obtained  from  the  open  literature. 

The  ab  initio  calculations  are  performed  using  the  CASTEP 
software  [60].  The  ultrasoft  pseudopotentials  are  used  in  conjunc¬ 
tion  with  the  Perdew-Burke-Ernzerhof  (PBE)  generalized  gradient 
approximations  (GGA)  exchange  correlation  function.  The  cutoff  of 
plane-wave  basis  set  is  500  eV  atom-1 .  The  energy  tolerance  for  the 
self-consistent  field  convergence  is  5.0  x  10-7  eV  atom-1  for  all  the 
calculations. 

The  optimized  parameters  obtained  from  the  procedures 
described  above  are  listed  in  Table  1.  Making  use  of  these 
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Table  2 

Physical  properties  of  5  LixSi  crystal  structures  using  the  new  2NN  MEAM  potentials.  The  units  of  the  cohesive  energy  (Ec),  equilibrium  lattice  constants  (A,  B  and  C)  and  elastic 
constants  (Cn-C66)  are  in  eV,  A  and  GPa,  respectively.  The  first  row  of  each  alloy  system  is  from  our  MD  simulations,  and  the  second  row  is  from  our  ab  initio  calculations. 


Ec 

A 

B 

C 

Cn 

c22 

C33 

C12 

Cl3 

C23 

C44 

C55 

C66 

Lii  Sii 

3.333 

9.346 

9.346 

5.673 

112 

112 

94 

58 

72 

72 

42 

42 

26 

3.301 

9.356 

9.356 

5.742 

101 

101 

77 

19 

37 

37 

45 

45 

37 

Li  1 5  Si4 

2.361 

10.589 

10.589 

10.589 

48 

48 

48 

33 

33 

33 

23 

23 

23 

2.466 

10.616 

10.616 

10.616 

47 

47 

47 

21 

21 

21 

28 

28 

28 

Li22Si5 

2.270 

18.776 

18.776 

18.776 

52 

52 

52 

28 

28 

28 

30 

30 

30 

2.353 

18.680 

18.680 

18.680 

46 

46 

46 

23 

23 

23 

35 

35 

35 

Li^Siy 

2.855 

8.139 

20.358 

14.343 

87 

78 

85 

45 

42 

39 

17 

20 

15 

2.950 

8.553 

19.647 

14.317 

92 

97 

90 

5 

11 

8 

28 

26 

24 

Lii3Si4 

2.412 

7.993 

15.103 

4.427 

75 

70 

74 

22 

24 

28 

20 

14 

14 

2.556 

7.932 

15.105 

4.437 

74 

61 

77 

17 

11 

10 

23 

24 

28 

parameters  in  Eq.  (1),  one  obtains  the  2NN  MEAM  interatomic 
potential  for  Li-Si  alloys. 

2.3.  Validity  and  accuracy  of  the  new  2NN  MEAM  potential  for 
LixSi 

To  assess  the  accuracy  of  the  new  potential,  we  have  calculated 
several  physical  properties  of  LixSi  alloys.  Shown  in  Table  2  are  the 
properties  of  crystalline  LixSi  alloys  calculated  by  ab  initio  simula¬ 
tions,  as  well  as  the  MD  results  by  using  the  2NN  MEAM  potential 
for  Li-Si.  It  is  seen  that,  using  the  2NN  MEAM  potential  in  con¬ 
junction  with  the  parameters  listed  in  Table  1,  the  MD  simulations 
predict  these  physical  properties  very  well. 

Another  comparison  is  shown  in  Fig.  1,  where  the  MD  results 
based  on  the  new  2NN  MEAM  potential  are  shown  together  with 
the  ab  initio  results  for  the  binding  energies  of  5  auxiliary  structures. 
Overall,  the  average  difference  is  estimated  to  be  less  than  10%, 
expect  for  the  “Dimer”  structure  which  gives  a  relatively  higher 
value,  ~60%. 

Another  important  validation  of  interatomic  potentials  for  crys¬ 
talline  materials  is  their  ability  to  simulate  the  transition  from  a  dis¬ 
ordered  structure  to  an  ordered  structure[52].  To  this  end,  MD  sim¬ 
ulations  are  conducted  using  the  new  2NN  MEAM  potential  to  simu¬ 
late  the  transition  of  several  Li-Si  alloys  from  disordered  structures 
to  crystalline  structures.  The  disordered  structures  are  obtained  by 
creating  an  offset  (ca.  0.75  A)  from  the  original  sites  of  the  atoms. 
Fig.  2a  shows  a  typical  disordered  structure  of  L^Sis.  Next  to  it,  the 
corresponding  radial  distribution  function  (RDF)  is  plotted,  which 
clearly  indicates  a  disordered  structure.  From  this  initially  disor¬ 
dered  structure,  MD  simulations  are  conducted  to  equilibrate  the 


Li-Si  nearest  neighbor  distance  (A) 


system.  During  the  equilibration,  the  structure  gradually  becomes 
more  and  more  ordered  and  eventually  returns  back  to  its  origi¬ 
nal  crystalline  structure  with  cohesive  energy  and  lattice  constants 
identical  to  the  original  ordered  structure.  Fig.  2b-d  shows  several 
snapshots  of  the  structure  and  the  corresponding  RDF  during  the 
equilibration.  Such  transition  can  be  considered  as  a  specific  tran¬ 
sition  path  from  the  molten  state  to  a  solid  state,  which  in  turn 
proves  that  the  newly  developed  2NN  MEAM  potential  is  robust. 

3.  Properties  of  amorphous  LixSi  alloys 

In  this  section,  properties  of  LixSi  alloys  will  be  computed  using 
MD  simulations  based  on  the  2NN  MEAM  interatomic  potential 
developed  in  the  previous  section.  The  first  step  in  MD  simulation 
is  to  construct  a  simulation  cell  with  atomistic  structure  identical 
to  that  of  the  material  being  simulated.  This  is  usually  easy  for  crys¬ 
talline  materials.  However,  it  is  not  so  trivial  for  amorphous  mate¬ 
rials.  In  this  work,  we  create  the  amorphous  structure  via  a  rapid 
quench  process  -  increasing  the  temperature  of  an  initially  crys¬ 
talline  structure  of  a  given  composition  to  4000  K;  then  decreasing 
the  temperature  rapidly  to  room  temperature.  Such  rapid  quench¬ 
ing  creates  an  amorphous  structure  with  desired  composition.  The 
amorphous  nature  of  the  resulted  simulation  cells  is  demonstrated 
by  the  RDF  shown  in  Fig.  3  for  all  LixSi  alloys  studied  in  this  work.  It  is 
clearly  seen  that  the  simulations  cells  created  by  the  rapid  quench¬ 
ing  are  indeed  amorphous.  In  all  these  alloys,  the  first  neighbor 
distance  varies  from  2.75  to  2.90  A  for  Li— Li  with  the  increment  of 
Li  concentration,  while  it  remains  constant  (z.e.,  2.67  A)  for  Li-Si. 
One  interesting  observation  is  that  the  Si  atoms  seem  to  be  much 
less  disordered  than  the  Li  atoms,  which  indicates  some  Si  atoms 
are  still  covalently  bonded  whereas  with  slightly  larger  bond  dis¬ 
tance  varying  from  2.35  to  2.54  A  after  lithiation.  Overall,  our  results 
are  in  good  agreement  with  the  recent  ab  initio  studies  [51,61  ]. 

Once  the  simulation  cell  is  constructed,  MD  simulations  can  be 
performed.  Unless  stated  otherwise,  MD  simulations  in  this  paper 
are  conducted  on  the  NPT  ensemble  at  finite  temperature  by  using 
the  Nose-Hoover  thermostat  and  Parrinello-Rahman  pressostat. 
The  cutoff  distance  is  set  up  as  4.8  A.  Time  step  is  set  as  1  fs.  Each  sim¬ 
ulation  starts  with  lOOOps  equilibrium,  followed  by  an  additional 
2000  ps  for  data  collection. 

3.  t.  Coefficient  of  compositional  expansion  (CCE)  of  amorphous 
LixSi  alloys 

The  coefficient  of  compositional  expansion  (CCE)  due  to  Li  inser¬ 
tion  can  be  calculated  by  the  linear  strain  per  lithium  concentration 
[62,63], 


Fig.  1.  Comparison  between  our  MD  simulations  (data  points)  and  ab  initio  calcula¬ 
tions  (solid  curve).  Energies  of  several  crystals  and  molecular  structures  are  plotted 
versus  the  Li-Si  nearest-neighbor  distance. 


x=x0 


(8) 
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Li-Si  distance  r  (A)  L'-S'  distance  r  (A) 


Fig.  2.  Several  typical  snapshots  of  intermediate  structures  and  corresponding  RDFs  during  the  transition  from  disordered  (a)  to  ordered  states  (d). 


where  Si  is  the  linear  compositional  expansion  and  x  is  the 
lithium  concentration.  The  strain  of  an  amorphous  structure  aris¬ 
ing  from  the  solute  concentration  change  is  purely  volumetric  and 
isotropic.  Thus  the  linear  strain  of  amorphous  Li*Si  alloy  can  be 
written  as 


vw-no) 

3V(0) 


(9) 


where  Vo  is  the  volume  of  crystalline  silicon.  In  this  paper,  the 
volume  change  V(x)  is  computed  by  conducting  MD  simulations 
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using  the  2NN  MEAM  potential  developed  in  the  previous  section. 
Shown  in  Fig.  4  is  the  computed  volume  expansion  versus  lithium 
concentration  in  LixSi.  The  results  agree  well  with  the  ab  initio  calcu¬ 
lations  and  experimental  measurements.  For  example,  for  x  =  3.75, 
our  result  is  3.9,  while  the  ab  initio  value  is  4.1  [42],  and  the  experi¬ 
mental  data  ranges  from  3.6  to  5.4  [41  ].  Making  use  of  these  results 
in  Eq.  (8)  yields  r)  =  0.2555. 


Fig.  4.  Volume  change  versus  lithium  concentration.  Experimental  results  are  eval¬ 
uated  from  [41]. 


3.2.  Elastic  modulus  of  amorphous  LixSi  alloys 

Elastic  constants  of  amorphous  LixSi  alloys  are  calculated  from 
MD  simulations  by  using  the  stress  and  strain  fluctuation  formula 
[64].  The  results  show  that  the  calculated  elastic  stiffness  ten¬ 
sor  is  isotropic,  as  expected  from  the  amorphous  structure.  Four 
elastic  constants  (only  two  of  them  are  independent)  are  shown 
in  Fig.  5.  For  pure  amorphous  Si  (x  =  0)  at  room  temperature,  the 
elastic  modulus  of  our  MD  result  is  116GPa,  slightly  larger  than 
the  experimental  measurement  (90-1 10  GPa)  [65],  while  the  shear 
modulus  (46  GPa)  is  within  the  range  of  available  experimental 
values  (45-56  GPa)  [65].  For  the  amorphous  alloys,  neither  exper¬ 
imental  nor  simulated  room  temperature  data  are  available  in  the 
literature.  For  comparison,  we  also  calculate  the  modulus  at  OK, 
which  is  plotted  in  Fig.  5  (solid  black  squares)  against  the  cor¬ 
responding  ab  initio  results  (open  diamonds)  from  [20].  It  seems 
that  our  MD  simulations  predict  stiffer  bulk  and  Young’s  moduli 
than  the  ab  initio  calculations  from  [20],  especially  for  pure  amor¬ 
phous  Si.  We  do  note  here,  however,  that  our  OK  modulus  for  pure 
amorphous  Si  agrees  well  with  other  ab  initio  studies,  e.g.  [66]  and 
empirical  calculations  [67]. 

To  use  the  concentration  dependent  modulus  in  a  continuum 
thermodynamic  model,  the  following  expressions  for  the  Young’s 
modulus  E(x)  and  the  Poisson’s  ratio  v(x)  may  be  used, 

E(x)  =  E(0)(1  +  r]EiX+r]E2X2),  v(0)  =  v(0)(l  +  r)v-ix  +  r]v2x2), 

(10) 

where  F(0)  =  102.6  GPa  is  the  Young’s  modulus  of  pure  amor¬ 
phous  Si,  t]ei  =  ~  0-41 ,  rjE2  =  0.047,  v(0)  =  0.30  is  the  Poisson’s  ratio  of 
pure  amorphous  Si,  r)v\  =0.16,  and  rjV2  =  -  0.018.  These  values  are 
obtained  by  fitting  the  quadratic  equations  to  our  MD  data  shown 
in  Fig.  5. 

3.3.  Composition-dependent  diffusivity  of  Li  in  LixSi  alloys 

The  diffusivity  of  Li  in  the  LixSi  alloys  depends  on  the  com¬ 
position  of  the  alloys  [68].  To  investigate  such  dependence,  MD 
simulations  are  conducted  to  calculate  the  mean  square  displace¬ 
ment  (MSD).  MSD  has  been  extensively  used  to  describe  the 
movement  of  atoms  in  solids,  liquids  and  gases.  Obviously,  the  MSD 
also  contains  information  on  the  diffusion  of  atoms.  It  has  been 
shown  that  the  MSD  increases  linearly  with  time  when  diffusion 
occurs  [69],  and  the  slop  of  the  MSD  versus  time  curve  gives  the 
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Fig.  5.  Elastic  constants  of  five  LixSi  amorphous  structures  obtained  from  the  new  2NN  MEAM  potentials  at  0  K  and  300 1<.  Dotted  lines  are  used  for  visual  convenience. 
Experimental  results  are  from  [65]. 


diffusion  constant  D, 

D  =  -  iim  Msp(^  +  A0-MSD(0 
6  At^oo  At 

1  /  |r,-(t  +  At)  -  r,(t)|2\ 

=  -  lim  - - - - — , 

6 1 — >oo  At 


(11) 


where  r,(t)  is  the  position  vector  of  the  diffusing  atom  i  at  the  time 
t,  and  the  symbol  ()  denotes  the  average  over  all  the  diffusing  atoms 
in  the  simulation  cell.  Theoretically,  the  MSD(t)  is  a  linear  function 
of  t  so  that  the  D  calculated  from  Eq.  (11)  is  independent  of  t.  In 
reality,  there  are  always  fluctuations  in  the  MSD(t)  when  calculated 
by  MD  simulations.  Thus,  linear  regression,  instead  of  numerical 
derivative,  is  used  to  calculate  the  diffusivity  D  from  the  MSD(t) 
curve. 

Another  challenge  is  computing  MSD  at  room  temperature. 
Because  the  diffusion  of  Li  at  room  temperature  is  rather  slow,  it 
is  extremely  time  consuming  to  conduct  MD  in  order  to  calculate 
MSD,  particularly  at  lower  Li  concentrations.  Therefore,  we  com¬ 
pute  the  MSD  at  elevated  temperature  in  this  work.  An  example 
is  shown  in  Fig.  6  where  the  MSD  of  Li  atoms  in  L^Sii  is  plot¬ 
ted  as  a  function  of  t  at  several  elevated  temperatures.  From  these 
high  temperature  data,  diffusivity  at  temperatures  can  be  obtained 
from  (1 1 ).  Once  the  diffusivity  is  known  at  several  higher  tempera¬ 
tures,  its  room  temperature  value  can  be  extrapolated  through  the 
following  Arrhenius  equation, 

D(T)  =  D  oexp(^),  (12) 

where  T is  the  absolute  temperature,  R  is  the  universal  gas  constant, 
and  D0  and  Ea  are  constants  independent  of  temperature.  In  fact, 


the  high  temperature  data  show  a  linear  relationship  between  ln(D) 
and  (1/T),  as  shown  by  the  solid  dots  in  Fig.  7. 

Using  the  method  described  above,  we  have  calculated  the  room 
temperature  (300  K)  diffusivity  of  Li  and  Si  in  both  crystalline  and 
amorphous  LixSi  for  several  values  of  x.  The  results  are  shown  in 
Fig.  8a  for  Li  and  Fig.  8b  for  Si.  Several  observations  can  be  made. 
First,  the  diffusion  of  Li  in  LixSi  alloy  is  at  least  an  order  of  magnitude 
faster  than  that  of  Si,  indicating  that  Li  is  the  dominant  diffusing 
species.  Second,  the  diffusivity  of  Li  either  or  Si  is  orders  of  mag¬ 
nitude  smaller  in  crystalline  LixSi  than  that  in  the  corresponding 
amorphous  alloy,  with  the  exception  of  Li  in  Lil2Si7,  where  the 
difference  is  much  smaller.  This  is  supports  the  hypothesis  used 
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Fig.  6.  MSD  of  Li  atoms  in  a  LiiSii  simulation  cell  under  different  temperatures. 
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Fig.  7.  Diffusivity  is  related  to  temperature  through  the  Arrhenius  equation.  The 
solids  dots  are  calculated  directly  from  the  high  temperature  MSD  data. 

in  Ref.  [70]  that  the  large  difference  in  diffusivity  between  the 
amorphous  and  the  crystalline  structures  is  a  major  cause  of  the 
observed  radial  cracking  of  Si  particles  and  wires  [39].  Third,  in 
the  amorphous  LixSi  alloys,  diffusivity  of  Li  increases  dramatically 
and  monotonically  as  the  Li  concentration  increases;  in  the  crys¬ 
talline  LixSi  alloys,  however,  the  diffusivity  of  Li  first  increases  with 
increasing  Li  concentration,  then  starts  to  decrease  as  Li  concentra¬ 
tion  increases  further.  The  highest  diffusivity  of  Li  seems  to  occur 
in  Lii3Si4.  Finally,  we  note  that  similar  trend  can  also  be  observed 
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Fig.  8.  Diffusivity  of  (a)  Li  and  (b)  Si  atoms  in  both  amorphous  and  crystalline  LixSi 
alloys  at  300  K. 


Fig.  9.  Uniaxial  stress-strain  curves  of  amorphous  Li*Si  alloys  at  300  K  under  the 
strain  rate  of  108. 

for  the  diffusivity  of  Si,  although  the  highest  diffusivity  of  Si  occurs 
in  Li12Si7. 

It  is  noted  that  the  Li  diffusivity  in  amorphous  LiiSii  shown  in 
Fig.  8  is  smaller  than  the  recent  ab  initio  calculation  [61]  in  which 
the  mixing  process  is  studied.  In  addition,  our  Li  diffusivity  is  much 
higher  than  the  experimental  data  [68],  and  shows  a  monotonic 
increase  with  increasing  x,  while  the  experimental  data  seem  to 
exhibit  a  “W”  shape  with  two  minimum  regions  at  Li2.i±o.2Si  and 
Li3.2iO.2Si.  These  differences  might  be  attributed  to  the  possible 
coexistence  of  amorphous  and  crystalline  microstructures  in  the 
partially  lithiated  experimental  samples.  Further  studies  are  being 
conducted  to  reconcile  the  differences. 

3.4.  Uniaxial  stress-strain  relationship  at  300 1< 

To  understand  the  mechanical  behavior,  MD  simulations  at 
300  K  are  conducted  to  obtain  the  room  temperature  uniaxial 
stress-strain  relationship  for  amorphous  LixSi  alloys.  The  simula¬ 
tions  using  the  NPT  ensemble  are  carried  out  on  an  amorphous 
LixSi  supercell  by  applying  a  constant  strain  rate  uniaxially.  Periodic 
boundary  conditions  (PBC)  are  prescribed  in  all  three  directions  of 
the  simulation  cell.  Traction  on  the  surfaces  normal  to  the  loading 
direction  is  set  to  zero  to  ensure  that  the  deformation  is  uniaxial. 
To  explore  the  strain  rate  effect,  three  rates,  107, 108,  and  109,  are 
used. 
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Table  Al 

Optimized  parameters  in  the  2NN  MEAM  potential  for  Li  [52]  and  Si.  The  units  of  the  cohesive  energy  Ec  and  equilibrium  nearest  neighbor  distance  re  are  in  eV  and  A, 
respectively. 


Ec 

re 

a 

d 

A 

£(°) 

pu 

£(2) 

£<3) 

f(2) 

t(3) 

Cmax 

Cmin 

Li[52] 

1.65 

2.99 

3.00 

0.14 

0.64 

1.03 

4.88 

4.15 

5.27 

-1.46 

4.13 

-0.57 

1.91 

0.31 

Si 

4.63 

2.35 

4.90 

0.03 

0.53 

3.00 

7.50 

0.00 

3.00 

1.45 

7.61 

-2.10 

2.60 

0.75 

Si[45] 

4.63 

2.35 

4.92 

0.00 

0.58 

3.55 

2.50 

0.00 

7.50 

1.80 

5.25 

-2.61 

2.80 

1.41 

Shown  in  Fig.  9  are  the  simulated  room  temperature  (300 1<) 
stress-strain  curves  under  the  strain  rate  of  108  for  various  compo¬ 
sitions  of  amorphous  Li*Si.  It  is  seen  that  the  stress-strain  curves 
exhibit  the  typical  characteristics  of  metallic  materials  -  an  ini¬ 
tial  linear  portion,  then  nonlinear  portion.  The  Young’s  modulus 
can  be  obtained  from  the  slop  of  the  linear  portion  of  the  uniaxial 
stress-strain  curve.  The  Young’s  modulus  so  obtained  agrees  well 
with  that  given  in  Fig.  5,  which  is  obtained  by  using  the  stress  and 
strain  fluctuation  formula  [64]  via  MD  simulation. 

The  yield  strength  can  then  be  extracted  from  the  uniaxial 
stress-strain  curves  in  Fig.  9.  This  is  done  by  first  plotting  the 
linear  regression  of  the  linear  portion  of  the  stress-strain  curve; 
then  plotting  a  parallel  line  with  0.2%  strain  offset.  The  intersection 
between  this  0.2%  offset  line  and  the  stress-strain  curve  gives  the 
yield  strength.  The  yield  strength  so  obtain  is  shown  in  Fig.  10  as  a 
function  of  Li  concentration  under  three  different  strain  rates. 

Clearly,  the  yield  strength  is  rate  dependent.  However,  the  data 
seem  to  show  that,  even  with  2  orders  of  magnitude  difference  in 
strain  rate,  the  yield  strength  does  not  change  significantly,  albeit 
this  observation  is  made  based  on  data  at  extremely  high  strain 
rates.  To  obtain  accurate  quasi-static  yield  strength,  more  studies 
are  needed. 

4.  Summary  and  concluding  remarks 

A  2NN  MEAM  interatomic  potential  for  the  Li-Si  system  is  devel¬ 
oped  in  this  paper  by  fitting  the  13  adjustable  parameters  in  the 
potential  function  to  a  number  of  physical  properties  of  the  Li-Si 
alloys.  These  properties  are  obtained  by  using  the  first  principles 
calculations,  and  the  fitting  is  done  by  using  the  particle  swarm 
optimization  method. 

Validity  and  accuracy  of  the  new  2NN  MEAM  potential  is 
demonstrated  by  computing  primary  structural  properties  for  both 
crystalline  and  amorphous  Li*Si  alloys,  as  well  as  simulating  the 
transition  from  disordered  to  ordered  states  of  the  atomistic  struc¬ 
ture. 

The  validated  2NN  MEAM  potential  is  then  used  to  calculate  sev¬ 
eral  thermomechanical  properties  of  the  Li-Si  systems  including 
the  elastic  modulus,  diffusivity  of  Li  in  Li*Si  alloys,  and  the  plas¬ 
tic  yield  strength.  The  results  show  that  these  properties  are  all 
concentration  dependent,  i.e.,  they  are  function  of  x  in  LixSi. 

This  new  interatomic  potential  will  be  powerful  tool  for  the 
modeling  and  simulation  of  fracture  failure  of  Si-based  anode  in 
Li-ion  batteries.  It  enables  us  to  validate  the  various  continuum 
models  on  Li  insertion;  it  can  be  used  to  estimate  some  of  the  mate¬ 
rial  properties  that  are  difficult  to  obtain  experimentally;  and  it 
gives  us  a  mean  to  study  the  microstructural  mechanisms  of  frac¬ 
ture  and  failure. 
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Appendix  A.  2NN  MEAM  potential  for  Si 

Using  the  same  method  described  in  Ref.  [52],  we  have  devel¬ 
oped  a  2NN  MEAM  potential  for  pure  Si.  The  parameters  in  this 


Fig.  Al.  Comparison  between  the  MD  results  (data  points)  and  ab  initio  calculation 
(solid  curve).  Energies  of  several  crystals  and  molecular  structures  are  plotted  versus 
the  Si-Si  nearest-neighbor  distance. 

2NN  MEAM  potential  are  somewhat  different  from  those  given  in 
Ref.  [45],  see  Table  Al.  The  elastic  constants  and  vacancy  energies 
calculated  by  our  2NN  MEAM  potential  are  consistent  with  the  val¬ 
ues  given  by  Baskes  [43],  Lee  [45]  and  Timonova  [71].  However, 
our  MEAM  potential  provides  a  better  description  of  the  relation 
between  energy  and  Si-Si  nearest  neighbor  distance  for  more  refer¬ 
ence  structures,  as  illustrated  in  Fig.  Al,  where  the  total  energies  of 
the  auxiliary  structures  are  plotted  versus  nearest  neighbor  distance 
varying  from  2.0  to  3.5  A. 
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